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Abstract. - The transport of magnetization is analyzed for the classical Heisenberg chain at and 
especially above the isotropic point. To this end, the Hamiltonian equations of motion are solved 
numerically for initial states realizing harmonic-like magnetization profiles of small amplitude and 
with random phases. Above the isotropic point, the resulting dynamics is observed to be diffusive 
in a hydrodynamic regime starting at comparatively small times and wave lengths. In particular, 
hydrodynamic regime and diffusion constant are both found to be in quantitative agreement with 
close-to-equilibrium results from quantum 5=1/2 autocorrelations at high temperatures. At the 
isotropic point, the resulting dynamics turns out to be non-diffusive at the considered times and 
wave lengths. 
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Transport in one-dimensional systems has been a topic 
of theoretical investigations for several decades, and even 
owadays there is an ongoing and still increasing interest 
in understanding transport phenomena in such systems, 
including their dependence on temperature [Tr[2"T]. Much 
work has been devoted to a qualitative classification of 
the dynamics into either non-normal ballistic or normal 
diffusive behavior. In this context the crucial mechanisms 
for the emergence of pure diffusion have been addressed 
and non-intcgrability is frequently discussed w.r.t. its role 
as a basic prerequisite, see jl], for instance. Particularly, 
spin chain models are a central issue of research [iH2"7]. 
with a considerable focus on transport of magnetization 
(spin) in the anisotropic S = 1/2 Heisenberg chain [IH2"2"]. 
which serves as a prototype model due to its integrability 
in terms of the Bethe ansatz [55] ■ The specific question of 
spin transport in this model has experienced an upsurge 
of interest during the last decade, not least motivated by 
experiments on low-dimensional quantum magnets, where 
large magnetic heat conduction [2l)H5T] and long nuclear 
magnetic relaxation times [551155] have been observed. But 
still, despite much effort, spin transport in the anisotropic 
5=1/2 Heisenberg chain is not fully understood yet, not 
even for finite anisotropies A/0. 

So far, many works have focused on the dynamics of 
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magnetization in the limit of both zero momentum and 
frequency in connection with the spin Drude weight, see, 
e.g., |5HS]. While non-zero Drude weights for A < 1 [8] im- 
ply partially ballistic dynamics [3] (for A = completely 
ballistic dynamics) and prevent the occurrence of purely 
diffusive dynamics at zero momentum and frequency, they 
do not allow for conclusions on diffusion laws at finite mo- 
mentum. Early analysis of the time-dependent correla- 
tion function of the local spin density [TU] suggested the 
absence of diffusion for all anisotropies < A < 1 in 
the high-temperature limit T = oo, see also |llj ; how- 
ever, recent low-temperature studies at the isotropic point 
A = 1, using bosonization and transfer matrix renormal- 
ization group [3] as well as Monte Carlo [T5], are consis- 
tent with finite-frequency diffusion at small momentum. 
On the other hand, T = oo investigations at A = 1 [15] 
on the basis of the Lindblad quantum master equation 
point to superdiffusive dynamics. But diffusion is found in 
this approach for anisotropics A > 1 [T3UT5] with a high- 
temperature diffusion coefficient which is in quantitative 
agreement with time-dependent perturbation theory for 
the current autocorrelation in the limit of zero momen- 
tum [TMTT] . Further analysis of the spin density autocor- 
relation at A > 1 [18] demonstrated finite-time diffusion 
at small momentum also for lower temperatures T < oo; 
however, the spectrum of the zero-momentum current au- 
tocorrelation might still feature low-frequency anomalies 
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indicative for an insulator [T51l20j . 

In this work we will shed light on the high-temperature 
dynamics of magnetization in a different way, namely, by 
applying classical mechanics. Certainly, the application of 
classical mechanics is not a novel concept; however, most 
work has focused on correlation functions at the isotropic 
point A = 1 yet 23 27 without providing a quantitative 
comparison to quantum mechanics and S = 1/2, which 
cannot be considered as a classical limit a priori, even at 
high temperatures T = oo. Here, we will provide such 
a comparison, motivated by a possibly simple scaling of 
correlation functions with S |16j . and we will particularly 
unveil a convincing quantitative agreement at A > 1. Due 
to this agreement, our results suggest that signatures of 
quantum 5=1/2 diffusion at A > 1 [TBTfTB] can be indeed 
understood "classical" phenomenon. 

The present paper is structured as follows: First of all, 
we will give the pertinent definitions by introducing the 
anisotropic Heisenberg chain in general and in particular 
the corresponding Hamiltonian equations of motion in the 
case of classical mechanics. Thereafter we will establish a 
transport scenario by specifying a certain class of initial 
states with inhomogeneous magnetization profiles of small 
amplitude. We will afterwards continue by discussing the 
properties of diffusive dynamics and in this context we will 
briefly review close-to-equilibrium results on diffusion at 
high temperatures in the case of quantum mechanics and 
S = 1/2. To these results we will eventually compare our 
numerical findings on classical mechanics, starting above 
but also proceeding to the isotropic point. We will finally 
close with a summary. 

In this work we investigate the anisotropic Heisenberg 
chain (XXZ model) with a Hamiltonian of the form 
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AS r *S* +1 ), (1) 



where S l r (i = x, y, z) arc the components of the spin 
S r at site r, L denotes the number of sites, J represents 
the exchange coupling constant, and A is the anisotropy 
parameter. 

In the case of classical mechanics (CM) the spins S r 
are three-dimensional vectors of unit length, |S r | = 1, and 
their dynamics is generated by the Hamiltonian equations 
of motion 



dH 
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with J r ±i = ± J S r x (£j? ±1 , Sf ±1 , A S^ii) as the incoming 
and outgoing spin currents. Since this set of equations is 
non-intcgrable in terms of the Liouville- Arnold theorem 
34,35], exact analytical solutions can only be given for 
mostly trivial initial configurations. We thus solve the set 
of equations numerically using a 4th order Runge-Kutta 
algorithm with a fixed time step of St J = 0.01; however, 
the presented results in the work at hand will not depend 



on this particular choice of the time step, cf. 
In order to realize a transport scenario we choose a certain 
class of initial states with non-homogenous magnetization 
profiles. Concretely, these initial states read 



S r (0) 



cos(a r ) \J\ — A 2 cos 2 (qr) a 2 



sin(a r ) y/l — A 2 cos' 
A cos(qr) a r 



{qr) a 2 
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where the first term of the z-component is a cosine-shaped 
magnetization profile at a single momentum q = 2nk/L 
and with the total amplitude A G [0, 1]. The second term 
of the z-component introduces local amplitudes a r G [0, 1] 
drawn at random from a uniform distribution. This term 
hence adds contributions at other momenta q' ^ q, but 
the main contribution is still at q. The remaining x- and 
y-components essentially contain local phases a r G [0, 2tt] 
drawn at random from a uniform distribution again. We 
emphasize that Eq. ([3]) defines a subclass of all possible 
initial states, which is particularly evident for small total 
amplitudes A <C 1. Later we will choose a small A = 1/4 
in order to model "close-to-equilibrium" states. Due to the 
random choice of local phases a r , these states correspond 
reasonably to high temperatures. In fact, also correlation 
functions at high temperatures arc obtained by sampling 
over a set of completely random configurations, see [26] . 
for instance. 

Our aim is to analyze the decay of the g-mode in time, 
i.e., the relaxation of S* = Y^,r cos(qr) s%. For long chains 
of length L = 18000, as considered throughout this work, 
the initial value Sg(0) ~ N A/4 does not depend on the 
concrete realization of the random numbers in Eq. ([3]) any 
further; however, the relaxation of S q still differs from one 
realization to another. We therefore analyze the ensemble 
average M q = l/N J2n S q over N = 1000 different initial 
states in order to obtain the "typical" dynamics. One may 
expect that each S q is close to the ensemble average M q in 
the thermodynamic limit L — > oo. In fact, this expectation 
is supported because our results for M q at large q can be 
reproduced for L < 18000 if N > 1000. In other words, 
the shorter the chain the more initial states are needed to 
reach the ensemble average. 

The dynamics may be called diffusive if M q satisfies a 
(lattice) diffusion equation of the form 



M q {t) = -9 D M q {t) , <f = 2[1 - cos(g)] 
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with a diffusion constant D, which is independent of both 
time t and momentum q. Of course, a diffusion equation 
is expected to be valid solely in a hydrodynamic regime 
starting at sufficiently long t and small q, e.g., in the limit 
of t — » oo and q — > 0. Hence, in order to investigate the 
full time- and momentum-dependence of the dynamics, we 
reformulate Eq. (U) and introduce a generalized diffusion 
coefficient, which given by [T8] 
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This quantity is in the focus of the present work. It is well 
defined whenever M q takes on a non-zero value, otherwise 
it diverges. Such a divergence naturally emerges after the 
relaxation of M q due to "fluctuations" around zero, e.g., 
as arising from numerical errors; however, we study the 
dynamics until the relaxation of M q here. In any case, a 
divergence of D q (t) is a clear indication of non-diffusive 
behavior. It is worth to mention that Eq. (|S|) is equivalent 
to writing M q in the form 



M q {t) 
M q (0) 



exp 



dt'D q (t'} 



(0) 



which leads to a strictly exponential decay for a diffusion 
constant D q {t) = const. 

Before presenting results for the case of CM, let us also 
discuss the case of quantum mechanics (QM) in detail, to 
which we are going to compare afterwards. In that case, 
all spins become operators, their time arguments have to 
be understood with respect to the Heisenberg picture, and 
initial configurations are specified by a density matrix p, 
which is possibly a pure state. In fact, pure states have 
been considered in numerical simulations |21| using the 
time-dependent density matrix renormalization group, yet 
restricted to zero temperature. In principle the "typical" 
dynamics can be defined again by the ensemble average 
M q (t) = 1/N J2n T AS q {t)p}- Here, however, we decide 
to proceed in a different way and instead define M q as 
the density autocorrelation M q (t) = (S q (t)St q (0)) , where 
the angles denote the canonical equilibrium average at the 
inverse temperature (3 = I/T, set to /3 = 0. This way we 
are going to compare to linear response theory [35] at high 
temperatures, i.e., the quantity D q {t) in Eq. ([5]) becomes 
a generalized high-temperature diffusion coefficient close 
to equilibrium. 

We especially make use of the perturbation theory (PT) in 
[SHUT] for S = 1/2 and A > 1, which yields a prediction 
for D q {t) in the limit of q — > 0. Therein, D q ^ (t) is found 
to take on the constant value D q ^o J — 0.88/ A at a time 
scale t J > 3.0/ A. Moreover, this value has been brought 
into quantitative agreement with diffusion constants from 
complementary non-equilibrium bath scenarios [13H15] on 
the basis of the Lindblad quantum master equation. In 
addition, D q (t) « D q ^ (t) has been observed in [TH] at 
finite momenta q < 27r/9, at least at the treatable times 
t J< 10, using full exact diagonalization (ED) for a chain 
of length L = 18. In view of these findings, the starting 
point of the following CM simulation will be A > 1 and 
q = 2tt/9. 

Now we turn to our CM findings. In Fig. [T] we display 
the numerical result for the decay of M q for anisotropics 
A = 1.5, 2.0 at momenta q = 2n/9, it/9, as obtained for 
a small total amplitude A = 1/4. To repeat, we consider 
a chain of length L = 18000 and average over N = 1000 
initial configurations. Apparently, the decay of M q takes 
place at rather short time scales t J < 30. This decay is 
the faster the larger A and the larger q. For comparison, 
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Fig. 1: Relaxation of modes M q (t) at momenta q = 7r/9, 2n/9 
for anisotropies A = 1.5, 2.0 in a chain of length L = 18000, 
as obtained numerically from classical mechanics for the total 
amplitude A = 1/4 (non-solid curves). Numerical parameters: 
TV = 1000, St J — 0.01. For comparison, the decay of quantum 
S = 1/2 autocorrelations at p = is shown, as predicted by 
the q — > perturbation theory in [T7] (solid curves). 



we additionally indicate the relaxation of S = 1/2, (3 = 
density autocorrelations, using Eq. ©, D q (t) « D q ^o(t), 
as well as D q ^o(t) according to the PT in [T7] (depicted 
in Fig. (2). The obvious agreement in Fig. [T] is remarkable 
for several reasons: First, we compare the average decay 
of a density to the decay of a density autocorrelation. Sec- 
ond, this comparison is done between CM and QM, but 
generally the quantum problem can only be considered 
reliably as classical for a sufficiently large spin quantum 
number 5, even in the limit of high temperatures. In any 
case, the agreement usually requires a proper scaling fac- 
tor taking into account the different length of a spin in 
CM and QM. The latter point deserves closer attention 
since we have not introduced such a scaling factor, i.e., 
the spins in Eq. ([2]) are still three-dimensional vectors of 
unit length. We emphasize that this observation is not an 
artefact of the applied numerical integrator because later 
we will reproduce existing CM simulations at A = 1. 

Next we intend to classify the type of the dynamics in 
Fig. [IJ To this end, the time-dependence of the generalized 
diffusion coefficient D q (t) is displayed in Fig. [3] at fixed 
momentum q = tt/9 for anisotropics A = 1.5, 2.0. For a 
small total amplitude A = 1/4 we again find a convincing 
agreement with D q ^o(t) according to the PT in [T7]. The 
CM curve D n /g(t) not only takes on an approximately 
constant value close to D q ^ J = 0.88/ A, but also reaches 
this value at roughly the same point in time, namely, at 
t J « 3.0/A. Clearly, the agreement is better for A = 1.5 
in Fig. dla) than for A = 2.0 in Fig. HJb), where the CM 
curve £> w / g (f) slightly decreases for t J > 10. In order to 
give insight into the origin of this decrease, we also show 
D-n-/g(t) for larger total amplitudes A > 1/4. Apparently, 
the decrease is the more pronounced the larger A. Thus, 
since the parameter A adjusts the total amplitude of the 
initial magnetization profile, we identify the decrease as a 
"not-close-to-equilibrium" effect. Even though not shown 
here, the decrease becomes also more pronounced as A 
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Fig. 3: Time-dependence of the generalized diffusion coefficient 
D^/ 36 (t) for anisotropics A = 1.5, 2.0 in a chain of length 
L = 18000, as obtained numerically from classical mechanics 
for the total amplitude A — 1/4 and different averages over 
N = 1000 and 10000 initial states (solid curves). Remaining 
numerical parameter: StJ — 0.01. As in Fig. [21 the quantum 
S = 1/2, /3 = prediction of the q —¥ perturbation theory in 
[T7] is indicated (dashed curves). 



Fig. 2: Time-dependence of the generalized diffusion coefficient 
D^/ S (t) for anisotropies (a) A = 1.5, (b) 2.0 in a chain of length 
L = 18000, as obtained numerically from classical mechanics 
for different total amplitudes A = 1/4, 1/2, and 1 (non-solid 
curves). Numerical parameters: N = 1000, St J — 0.01. The 
q — > prediction [T7] for quantum S = 1/2 autocorrelations at 
P = is indicated (thick solid curves). In (b) the numerical 
q = 7r/9 result for the classical autocorrelation at j3 = is 
shown for N = 10000 (thin solid curve). 

increases at fixed A, which is already visible in Figs. [2][a) 
and (b). However, for A = 2.0, the effect is comparatively 
small for A = 1/4 and may completely disappear in the 
limit of A —> 0. In our CM simulation we cannot consider 
the case of a very small total amplitude 4 < 1/4 since 
D^/ 9 (t) develops strong "oscillations" if A is decreased to 
much. Such oscillations are already observable for A = 1/4 
in Fig. [5] and essentially arise from a too small number of 
initial states, as we will demonstrate later in detail. This 
observation indicates that averaging is more important for 
"close-to-cquilibrium" states, at least for a chain of finite 
length. It is worth to mention that M q is less sensitive 
to averaging because the integration in Eq. ((6]) smoothes 
the oscillatory behavior. For completeness, we illustrate 
in Fig. [5](b) the same oscillatory behavior for the classical 
density-density correlation at /? = 0, evaluated according 
to [55] and sampled over TV = 10000 completely random 
initial states (10 times more than before). Despite the still 
strong oscillations, this correlation is consistent with the 
decay at A = 1/4, which particularly excludes a possibly 
"far-from-cquilibrium" scenario for this choice of A. 

So far, the CM results in Figs. [1] [2] can be reproduced 
for a chain of length L = 18 if the averaging is performed 
over N ^> 10000 initial configurations, which is feasible 
for a short chain. Hence, these results are certainly free 
of relevant finite-size effects. Next we proceed to smaller 
momenta q < 7r/9, which cannot be realized for L = 18 



any further. In Fig. [3] we display the generalized diffusion 
coefficient D q {t) at fixed momentum q = 7r/36, still for 
anisotropies A = 1.5, 2.0 and the small total amplitude 
A = 1/4. The similarity to the results in Fig. [5] is obvi- 
ous, even though a wider time window t J < 120 is plotted 
in Fig. [31 In this time window, the decrease of D n / 36 (t) 
appears to stagnate. On that account, without causing a 
large error, we may consider D v / 36 (t) as an approximately 
time-independent quantity, which is not identical but at 
least close to D q ^ from the PT in [17] . We have checked 
in addition that the CM curve for D q (t) does not change 
significantly at q = n/18 or 7r/72, although not shown ex- 
plicitly here. This independence of q, together with the 
approximative independence of t, is a clear signature of 
diffusion. This observation is one main result of the work 
at hand. It indicates the existence of CM diffusion in an 
extended region of time and momentum in quantitative 
agreement with QM predictions for S = 1/2 at (3 = 0. 
Of course, even smaller q and longer t may be evaluated; 
however, it is reasonable that the above t- and g-region 
lies already within the hydrodynamic regime and at some 
point it is hardly possible to exclude errors due to the 
applied numerical integrator. In this context we demon- 
strate finally that "oscillations" of D q (t) do not refer to 
such errors but only depend on the number of initial con- 
figurations. To this end, Fig. [3] also shows D^/ 36 (t) for 
N = 10000, which clearly is smoother than for N = 1000. 

Next we turn to the isotropic point, i.e., A = 1. First, 
we intend to clarify whether there still is a relationship to 
S = 1/2, (3 = density autocorrelations, as done before 
for A > 1 in Fig. \T\ However, since the PT in [T7] is not 
applicable in the region A ~ 1, we proceed differently and 
instead obtain these density autocorrelations numerically 
using ED for a chain of length L = 18. To repeat, in such 
a short chain the smallest realized (non-zero) momentum 
is q = 7r/9. In Fig. 0] the resulting relaxation of density 
autocorrelations is displayed for momenta q = 2-7r/9, 7r/9 
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Fig. 4: Relaxation of modes M q (t) at momenta q — it/9, 2n/9 
for the anisotropy A = 1.0 in a chain of length L — 18000, 
as obtained numerically from classical mechanics for the total 
amplitude A = 1/4 (non-solid curves). Numerical parameters: 
N — 1000, 8tJ = 0.01. For comparison, the decay of quantum 
S = 1/2 autocorrelations at /3 = is shown, as found from 
exact diagonalization for a chain of length L = 18. 



and compared with the CM result on the decay of M q at 
the same momenta but in a much longer chain of length 
L = 18000. Again, for a small total amplitude A = 1/4 
the quantitative agreement is remarkably good. Of course, 
the agreement appears to be much better for the larger 
momentum q = 2w/9, where the same slight oscillation 
around zero is found in CM and QM, emerging at a time 
scale tJ^Q. Such oscillations are well known to emerge 
at large q [55] and zero-crossings yield divergencies of the 
generalized diffusion coefficient [IS] , as also visible in later 
numerical data. For the smaller momentum q = ir/9 the 
CM and QM curves in Fig. [4] start to deviate from each 
other at a time scale t J ~ 10; however, finite-size effects 
of a QM curve for L — 18 usually set in at this time scale 
[TB] . On that account the agreement in Fig.H]supports the 
relationship between CM and QM, which may also hold 
in the limit of smaller momenta and longer times. 

In Fig. [5] we eventually summarize our CM results on 
the generalized diffusion coefficient D q (t) for A = 1 over a 
wide range of momenta q = 2 1 ~ l it/9 (i = 0,1,..., 5), i.e., 
between 2n/9 and 7f/144. (Strictly speaking, the smallest 
momentum reads q = 2ir 63/18000 « 7r/144.) Apparently, 
D q (t) depends significantly on momentum at q < 717 '72 
and even diverges because of the slight oscillation of M q 
around zero, as already shown in Fig. 2J In the depicted 
time window, t J < 120, a negligible dependence of D q (t) 
on momentum is found firstly when varying q from 7r/72 
to 7t/144. Although all CM curves in Fig. [5] arc consistent 
with D q (t) r-j D v im(f) + q A (t J) 2 /4, this (/-scaling has to 
be taken as a rule of thumb for the dominant dependence 
on q without derivation. We emphasize that the CM curve 
at q = 7r/144 agrees well with previous results on CM 
density-density correlations. For instance, using in Eq. JS]) 
the function exp[-0.537<? 2 (l + 0.1 ln|<?|) Jt ln(Jt)], which 
is one possible fit function in [5(5] (see also the discussion 
in [57]), we obtain for D n / 144(f) a corresponding function 
0.33 [1 + \n(t J)], as indicated in Fig. [5] (cirlces). Clearly, 
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t [units of 1/J] 

Fig. 5: The generalized diffusion coefficient D q (t) at momenta 
q — 2 1 ~ l tv/9 (i — 0, 1, . . . , 5) for the anisotropy A = 1 in a 
chain of length L — 18000, as found numerically from classical 
mechanics for the total amplitude A = 1/4 (using N = 1000 
initial states and a time step of St J = 0.01). For comparison, 
the function 0.33 [1 + ln(t J)] is indicated (circles). 



the time-dependence of D^/ 144 (t) may be a pointer to such 
a logarithm. Furthermore, because of the good agreement 
with [26], this logarithmic increase cannot be understood 
as a "not-close-to-equilibrium" effect resulting from the 
choice of the total amplitude A — l/A. In fact, we do not 
find a relevant dependence on A. One might be tempted 
to speculate on the time- and momentum-dependence of 
D q (t) beyond the depicted time window, see also a related 
discussion in |26827| : however, it is hardly feasible to verify 
a possibly logarithmic time-dependence of D q (t) from a 
merely numerical simulation. In any case, at least for the 
considered t and q in Fig. [5] the dynamics turns out to 
be non-diffusive. In view of the agreement with S = 1/2, 
P — density autocorrelations in Fig. [4] we may expect 
a similar dynamics in QM. This expectation has impact 
on numerical approaches to QM dynamics, simply since 
most approaches are confined to the here considered t and 
q, e.g., q = 7r/144 represents the smallest momentum in a 
long chain of length L = 288. While this length will never 
be reached by ED, it may be treatable by time-dependent 
density matrix renormalization group [21] . even though 
restricted to zero temperature yet. It is worth to mention 
that a recent result for L = 256 from a non-equilibrium 
bath scenario on the basis of the Lindblad quantum master 
equation supports supcrdiffusivc transport at ft = |13j . 
which is consistent with Fig. [5] 

In summary, we studied the transport of magnetization 
for the classical Heisenberg chain at and especially above 
the isotropic point. To this end, we solved the Hamiltonian 
equations of motion numerically for initial states realizing 
harmonic-like magnetization profiles of small amplitude 
and with random phases. Above the isotropic point, i.e., 
for large anisotropics A ~ 1.5 we found that the resulting 
dynamics of the ensemble average becomes diffusive in a 
hydrodynamic regime starting at rather small times and 
large momenta. We further brought the diffusion constant 
into good quantitative agreement with quantum 5 = 1/2 
predictions at high temperatures (3 = 0. At the isotropic 
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point, i.e., for the anisotropy A = 1 we did not observe 
evidence of diffusive dynamics at the considered times and 
momenta. But even for A = 1 we found at large momenta 
a remarkably good agreement with numerical results for 
quantum S — 1/2 density autocorrelations at /3 = 0. On 
that account we finally argued that our findings for the 
classical dynamics may also have impact on the quantum 
dynamics for A = 1 . While it is questionable if the present 
classical approach to transport is indeed appropriate to 
simulate possibly pure and strongly S'-dcpendcnt quan- 
tum effects at low temperatures, the possibility of a clas- 
sical simulation at not too low temperatures appears to 
be promising and may be explored further. 

* * * 

The author sincerely thanks W. Brenig, P. Prelovsek, 
J. Hcrbrych, T. Prosen, and M. Znidaric for interesting 
and fruitful discussions. 
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